Statistical evaporation of rotating clusters. II. Angular momentum distribution 
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The change in the angular momentum of an atomic cluster following evaporation is investigated 
using rigorous phase space theory and molecular dynamics simulations, with an aim at the possible 
rotational cooling and heating effects. Influences of the shape of the interaction potential, anhar- 
monicity of the vibrational density of states (DOS), and the initial distribution of excitation energies 
are systematically studied on the example of the Lennard- Jones cluster LJ14. For this system, the 
predictions of PST are in quantitative agreement with the results of the simulations, provided that 
the correct forms for the vibrational density of states and the interaction potential are used. The 
harmonic approximation to the DOS is used to obtain explicit forms for the angular momentum dis- 
tribution in larger clusters. These are seen to undergo preferential cooling when thermally excited, 
and preferential heating when subject to a strong vibrational excitation. 

I. INTRODUCTION 

When a polyatomic system is vibrationally excited, it may undergo spontaneous dissociation on a suitable time 
scale. Most often, different dissociation channels arc available, which may or may not be selected upon energetic 
criteria only. Rotational excitations are also very important, because they directly influence the relative amount 
of kinetic energy released in translational modes. Kinetic energies released are of constant use in experiments for 
e xtra cting information about binding energies, temperatures, and more generally nanocalorimetry data (see Refs. 
ITl U'il and references therein) . 

Rotational excitations practically act in two distinct ways. The first is direct, and reflects the greater ease of 
fragmenting the system due to centrifugal forces. The second way is indirect, and lies in the coupling between 
vibrational and rotational motion: because the system is not rigid, it may rearrange so as to reduce its angular velocity, 
thereby enhancing its stability. When an isolated cluster dissociates into two fragments, the angular momentum J is 
splitted into two parts, namely the orbital momentum L of the fragments in their relative motion, and the remaining 
angular momentum J r associated with the fragments themselves. It is very easy to see that, at low initial value of J, L 
nearly compensates J r . Hence, for a large system emitting a single atom, rotational heating is expected at low J. On 
the opposite, rapidly rotating clusters are more likely to lose angular momentum following evaporation, a phenomenon 
known as rotational cooling. These aspects of fragmentation have been previously investigated by Stace^ who used 
simple statistical models to emphasize the role of rotational temperature in the interpretation of experiments. 

In our previous work^ we performed molecular dynamics (MD) simulations as well as theoretical calculations based 
on the phase space theory (PST) to extract the distributions of kinetic energies released upon evaporation in small 
Lennard- Jones (LJ) atomic clusters. Our main result was to show that PST could be quantitatively accurate in a broad 
range of total energies and angular momenta, provided that the vibrational and rotational densities of states, as well 
as the interaction potential between the fragments, were all correctly accounted for. This result was previously proved 
by Weerasinghe and Amar for nonrotating clusters^ and partially observed for aluminium clusters by Peslherbe and 
Hase^*£ though these authors used a simple C jr n form for the dissociation potential, which induced some deviations 
between theory and simulation at large excitations. 

Here we will focus on angular momentum properties, namely their distributions and mean values, and how they are 
affected by the initial excitation. The present approach follows our previous effort^ with the same aim at assessing 
phase space theory in a quantitative comparison with simulations. In the next section, we briefly recall the main 
equations underlying the PST calculation of the angular momentum properties of the product cluster, and we treat 
separately the simple cases of an harmonic vibrational density of states, an interaction potential with the form C/r", 
or an initially nonrotating cluster. Our application to small and medium size Lennard- Jones clusters for various types 
of excitations is treated in Sec. IIIII In Sec. IIVI we discuss our results with respect to the works by Stace^ and those 
by Peslherbe and Hasej& before finally concluding. 



2 



II. PHASE SPACE THEORY 



In this Section, we give the main expressions for the angular momentum distribution of the product cluster following 
evaporation. To obtain such information, conservation of total angular momentum and energy has to be included, 
which is the case in the phase space theory developed by Pechukas and Lightj2iiS Klotsjii* 1 ^ and Chesnavich and 
Bowers. 13 i 14 i 15 

Here we consider a parent cluster characterized by a rotational angular momentum J and a total rovibrational energy 
E. We denote by J r the rotational angular momentum of the product cluster after dissociation. The unnormalized 
probability of finding a dissociation event characterized by the angular momentum J r of the product cluster within 
dJ r and the energy released £tr within de tr can be written asi^ 



V(e tI ,J r ;E,J) = Q^(E-E -e tT ) T(e* r ,J r )dL. (1) 

A mm (£„,J r ) 



In this equation, flif^ is the vibrational density of states of the product cluster having angular momentum J r , 
r(e*, J r ) is the available volume of states with rotational energy lower than e* in angular momentum space. For a 
given e tr , the rotational density of states is given by the sum of this quantity over the accessible area in the (J r , L) 
plane, with L the orbital angular momentum. This area is limited by mechanical (angular momentum and energetic) 
constraints, and will be discussed below. If both e tr and J r are fixed, the integration must be performed in a range of 
L limited by L m j n and L ma x- In this paper, we are interested in the distribution of angular momentum of the product 
cluster. The probability of finding a particular value of J r within dJ r is simply the integral of V over all possible 
values of e^ t : 

P{J r ;E,J)= [ V(£ t r, J r ;E,J)dStr (2) 

In PST the products are located at the transition state. Here we will consider rather large clusters, and our main 
assumption will be to treat the evaporative system LJ„ + i -^LJ„+LJ within the sphere+atom model. In this case, 
the available volume r(e*, J r ) is simply given by 2J r m& Therefore the probability of finding a dissociation event with 
angular momentum J r , starting from a total energy E and a total angular momentum J is given by: 

P{J r ;E, J) = P x 2J r / Slg*\E -E - e tr )[L max {e tI , J r ) - L min {e tT , J r )]ek tr , (3) 

where the normalization constant P accounts for channel and rotational degeneracies, and also depends on the parent 
density of states Vl^+i{E). In the following, we will make use of the notation AL(e tr , J r ) = £max(£tr, Jr) — -^min(£tr, Jr)- 
The problem is now reduced to determining A_L(e tr , J r )- 

We first consider the case of a radial dissociation potential V(r) given by V(r) = —C n /r p , with p > 2. Since the 
kinetic energy must be positive at the centrifugal barrier, the well known energetic constraint can be easily found: 

B n J 2 r +L 2 ^P-V/A p <e tI , (4) 

where we have introduced the rotational constant B n of the product cluster LJ„, and the Langevin parameter A p 
given as a function of p, C n , Ti and the reduced mass fi n by 

9 /// n \P/(P- 2 ) 

The second boundary comes from the conservation of angular momentum, J = J r + L, or more conveniently: 

\Jr - L\ < J < J r +L. (6) 

Let us denote by C the set of (J r , L) points that fulfill the equation e tr = L 2p ^ p ^ / A p + B n J 2 . For a given value of 
J r , we define £™ m (resp. ej."*) as the values of Etr, which correspond to the intersection between C and L = \ J — J r \ 
(resp. L = J + J r ). We find 

/_/ |2p/(p-2) (J+J ) 2 p/(p- 2 ) 
s^ = B n J 2 + l -^—!ll an d £ £t = Bn j2 + ki±^ (7) 
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Obviously, e™ m and e™* can not be larger than the maximum available energy E — Eq, and the largest value for J r will 
be reached when e^ ln =E — E . In Fig.^ a simple scheme summarizes the integration procedure in the (J r , L) plane. 
In this figure, C m i n and Ci nt are the contours C for e tr = e™ m and e tr = £tr\ respectively. When e"" n < e tr < e^, 
AL(e t r, J r ) is given by £(£tr, J r ) — \ J — J r \, where £(£tr> J r ) — [A p (etr — B n J 2 )](p~ 2 )/ 2 p. On the other hand, when 
eg* < etr < E - E Q , Ai(e tr , J r ) = (J + J r ) - \ J - J r \. 




FIG. 1: Schematic representation of the (J r ,L) integration plane. 



One interesting case is the harmonic limit of this model. For this we consider the harmonic VDOS for the product 
cluster, ft n (E) oc E s ~ 1 with s — 3n — 6. After replacing the expressions for AL(en, J r ) obtained above, equation © 
now becomes 



£ int 

P h (J r ;E,J) cx 2J r [ " (E - E a - e^y- 1 ^, J r )de tT 
- 2J r \J~J r \ J 11 (E - E - etrY^detr 



E-E 

s-1. 



+ 2J r (J + J r -\J- J r \) / {E-Eo- E tI ) s - L de tI , (8) 

By considering the analytical expression for £(et r , J r ) given previously and using the notation 7 = (p— 2)/2p, we thus 
obtain 

P h (J r ;E, J) oc 2J r A; f " (E-E - e tr ) s " 1 (et r - B n J^de tr 

- J ' " Jr kE-E -s tl .y]f:: 



s 

2J r (J + J r ~ \J-J r \) 

\E-Ea' 



(E-E -e ti y] E « (9) 



This equation can be further simplified into 



2 J r ( J + J r ) _int\s 



P h {J r ;E,J) cx (E — E — ) 5 

s 

_ 2J r [ J - J r | (Jg _ _ ggin)3 + 2 ^ (10) 
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where I sn is short for the integral in Eq. This integral can be written explicitely as 

i sn = g ( s - y;ZT^ E ~%-^ nJ?)k [ {j + j " )i+(s " )/7 - 1 j - j r\ i+{s - k)h ] ■ (id 

All the formalism above has been derived by assuming an interaction potential with the form V(r) = —C n /r p . As 
discussed in our previous work^ this expression does not give a very good account of the finite extension of the cluster, 
and a better representation of the atom-cluster interaction is provided by V(r) — —C n /(r — ro) p , with ro > 0. In 
this case, and more generally for an arbitrary form of V(r), the computation of £ m in(£tr, Jr) and £ m ax(etr, Jr) must 
be carried out numerically. For a series of L, the location r*(L) and the height £'(£) of the centrifugal barrier is 
obtained. At a given e tr , the integration boundaries become 



et(L) + B„J r 2 <6tr, 

\Jr - L\ < J < J r + L 



(12) 



and the limits L m in(etr, Jr) arL d L max (et r , Jr) can thus be directly calculated from the knowledge of £^(L). 

When we do not specify the harmonic limit to extract P(J r ; E, J), equation JSJl is used with the anharmonic 
vibrational densities of states f2( Jr )(i?), which includes some dependency with J r . To calculate fl^ Jr \E), we have 
used the Monte Carlo method proposed in Ref.0, further improved with the parallel tempering accelerating schemed 
The multiple histogram method^ was then used to estimate the configurational densities of states and, through a 
simple convolution product, the total VDOS. In practice, we have neglected the variations of £1 with J r , and used the 
approximation J r ~ J. As will be seen below, this does not entail a significant error. 

The dissociation potential V(r) felt by an atom leaving the n-atom LJ cluster was taken from our previous work^ 
and approximated by the form —C§/(r — ro) 6 form. The coefficients Cq and ro have been deduced from a fitting of 
constrained Monte Carlo simulations at finite temperature, as well as some Wang-Landau simulations. Their values 
are given in Ref. U 

For each value of angular momentum, we have also carried out one or several sets of 5000 independent molecular 
dynamics trajectories, used as reference data. The angular momentum distribution has been calculated for different 
values of the initial total angular momentum between J = and J ss 5 LJ units^ For the details of the MD 
simulations, we also refer the reader to Ref. 0. 



III. RESULTS 



The energetics of unimolecular evaporation has been considered in our previous article^ and we only focus here 
on the angular momentum distribution following the evaporation of a rotating atomic cluster. In the form detailed 
above, phase space theory directly gives us access to this distribution. The sharing between rotational energy of the 
product cluster and the translational kinetic energy released can then be analysed. 

We first consider the unimolecular dissociation involving the nearly spherical LJ13 product cluster. The reaction 
studied here is thus LJ14 — >LJi3+LJ. In Fig.[2we have plotted the angular momentum distribution of the LJ13 product 
cluster, as estimated from PST. The rotational constant B of LJ13 was taken at T — 0.3, such that B = 0.0444 LJ 
units£ Two values of the initial angular momentum, namely J — 0.5 and J = 3, were considered and the initial 
rovibrational energy was taken equal to E = —29 LJ units in both cases. Different approximations were assessed 
and a good overall agreement with MD results is obtained. In Fig. Ufa), two different interaction potentials were 
selected along with the anharmonic density of states for both calculations. The interaction potential corresponding 
to T — 0.3 is computed from MC simulations, and the ro = approximation refers to the long-range part of the 
exact LJ potential, with Cq — An units. As can be seen from Fig.[2Ja), the ro = approximation performs slightly 
less than the numerically exact radial potential, but still allows to reproduce quite accurately the MD results. The 
ro = approximation tends to underestimate the angular momentum of the product cluster, a feature previously 
noticed by Peslherbe and Hase^ In Fig.^b), the harmonic limit has been tested for J = 3. Even though a very good 
agreement is obtained when considering the anharmonic DOS, some discrepancy with MD becomes to appear in the 
harmonic case. A change in slope has to be noted at J r = J in Figs.J^a) and especially Hfb) . This discontinuity in 
the derivative of P(J r ; E, J) is due to the modulus \J r — J\ term in the definition of £™ m . Taking the derivative of 
Eq. |21 as a function of J r actually leads to the following value for the change in slope accross J = J r : 
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= -4P J / rt n J \E - E - e tI )de tI . (13) 
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FIG. 2: Angular momentum distribution of JV for the reaction LJ14 — >LJi3+LJ at E = —29, for two different values of the 
initial angular momentum J. (a) J = 0.5; (b) J = 3. 



Therefore the slope necessary decreases at J r = J, independently of the dissociation potential. In order to put this 
effect into evidence, we have run 15 sets of 5000 trajectories at J = 3.0 [see Fig.|3£b)]. Within this statistics, a change 
in the slope around J r — J can be effectively characterized, as seen in the statistical PST formalism. 
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FIG. 3: {J r } versus J for the reaction LJ14 — >LJi3+LJ at E — —29 within different approximations: (a) effect of the radial 
potential; (b) effect of the anharmonicity of the PES. The J — values, obtained from Eq. 1141 . are represented by stars or 
crosses. 



In Fig. |3| the mean angular momentum of the product cluster ( J r ) has been plotted as a function of J at the 
rovibrational energy E = —29. The same assumptions as in the previous figure have been tested. The line ( J r ) = J 
has also been shown in order to visualize more easily the rotational cooling or heating effects. First of all it is striking 
that the statistical PST model allows to perfectly reproduce the MD results when both the anharmonic DOS and the 
numerically exact radial potential are considered. The ro = approximation tends to underestimate (J r ) in the low 
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J regime. As J goes to 0, one recovers the exact limit of Eq. @ given by 

P{J r ;E, J = 0) = P x 2J r [ ° OW(£-£ - £tt )(fe tr . (14) 

The harmonic approximation to this equation is straightforward: 

P{J r ;E,J = 0) oc J r [E - E - B n J? - e f (J r )] S . (15) 

On the contrary, in the high J limit, (J r ) is overestimated. As J increases, the orbital transition state reaches 
small values, which explains that the ro = radial potential does not allow a correct description of the dissociation 
process. On the other hand the harmonic limit tends to underestimate (J r ) on the whole domain of J considered in 
this study. Therefore the harmonic model with ro = could well give results in quantitative agreement with MD in 
the low J regime, but this would remain fortuitous. 

Let J c be the value of the initial angular momentum for which (J r ) = J. For J < J c , evaporation leads to an 
increase in the mean angular momentum, the so-called rotational heating process. High initial values of J more likely 
result in a decrease of angular momentum corresponding to rotational cooling. At E — —29 LJ units, J c turns out 
to be approximately equal to 2.2 units. As seen from the previous figures, J c weakly depends on the dissociation 
potential, but is more sensitive to the harmonic approximation to the vibrational density of states. 
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FIG. 4: J c as a function of -E/n for the reaction LJ14 — >LJi3+LJ. 

In Fig. 0] we follow the evolution of J c as a function of the internal energy per atom E/n. J c was calculated with 
the numerically exact radial potential, both in the harmonic and anharmonic approximations to the VDOS. In the 
low energy regime, the cluster tends to be harmonic and both approximations give similar results. At higher energies, 
especially close to the melting phase change of LJ13 and beyond, stronger differences appear. Between E/n = 1.1 
and E/n — 1.4, J c reaches a plateau similar to the plateau in the microcanonical temperature at the same energy, 
but only within the anharmonic description. The same plateau must also be correlated with the one in the average 
kinetic energy released, as seen in Refs.(H and[l7l Hence, large internal energies in the parent cluster indeed requires 
an anharmonic treatment of vibration to extract quantitative informations about the threshold between rotational 
cooling and heating effects. 

Up to now we have only analyzed product angular momentum distributions for a given initial angular momentum J 
and an internal rovibrational energy E. We will now consider physically realistic conditions, where both internal energy 
and angular momentum are drawn from thermal distributions. The initial angular momentum distribution /(J) for the 
spherical cluster LJ„ + i is given by /(</) cx J 2 exp (— B n+ \ J 2 /^B^rot) • The distribution of the initial vibrational energy 
f(E) follows the VDOS, and we have allowed for a possible extra amount of energy AE brought by a specific excitation 
process (collision, photoexcitation, ionization,...). Therefore f(E) oc fl^^E — AE) exp [— (E — AE)/kBT v ^]. The 
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distribution of product angular momentum then reads 



P(J r ) oc U J 2 VL n+1 (E - A£)exp 



AE B n+1 J 2 



k R T, 



S J rot 



P(J r ;E,J)dJdE 



(16) 



The initial vibrational and rotational temperatures may be different, because the extra energy AE is converted into 
vibrational energy. In Fig. [5] two physical situations have thus been studied. The first one corresponds to a parent 





FIG. 5: Initial and final angular momenta following the evaporation of LJ14 for different physical conditions characterized by 
their thermal distributions and the extra excitation energy, (a) T = 0.05 and AE = 5; (b) T = 0.05 and AE — 8; (c) T = 0.20 
and AE = 0; (d) T = 0.40 and AE = 0. 



cluster prepared at low temperature (T = 0.05), for which two possible values of AE have been added, namely either 
AE = 4 or AE = 8 LJ units. This situation is typical of a cluster prepared in a supersonic expansion and photoexcited 
in the infrared domain. In Fig. Efa), the final distribution of J is slighty shifted to lower values. This is indicative 
of rotational cooling. In this case, the mean E/n is almost equal to 0.4. An extrapolation from Fig. 0] shows that 
J c will be approximately equal to 0.9, a value that is smaller than the initial angular momentum. It explains why 
cooling occurs under these physical conditions. On the other hand, in Fig. |S£b), E/n rs 0.7 and J c is increased to 
about 1.7. This value is larger than the mean value of the initial angular momentum, which is consistent with the 
observed rotational heating. 

We have also simulated a pure thermal excitation of the cluster at a temperature large enough to induce dissociation 
on a reasonable time scale, that is accessible to MD. Therefore the thermal distributions are directly considered for E 
and J with AE = 0, and T rot = T v ;b- We have considered two different temperatures, T = 0.2 and 0.4, and in both 
cases rotational cooling appears as the direct consequence of (J) > J c . The cfficcncy of the cooling process can be 
analysed by adjusting the product angular momentum distribution to a Boltzmann law. For example, the thermalized 
cluster initially at — 0.4 used in Fig.[SJd) yields a final rotational temperature close to T^~> — 0.16. At this initial 
temperature, E/n is approximately 1, and from Fig. 0] we see that the harmonic approximation would significantly 
underestimate rotational cooling. 

Finally we have investigated larger clusters, within the harmonic approximation as well as the simplest case for the 
radial potential, given by rg = and Cq = An LJ units. The rotational constant B n has been taken proportional 
to n -5 / 3 , with the proportionality factor calculated from the LJ13 cluster at T — 0. The value of E$, the energy 
difference between LJ n+ i and LJ ra , has been taken to 6.4 for all sizes, a good approximation in the present size range 
of interest. 

In Fig. the angular momentum distribution before (solid line) and after (dashed line) evaporation is shown for 
n = 50 and n = 100 with an initial (vibrational and rotational) temperature T = 0.5. The adjustment to a thermal 
distribution is also shown. At this initial temperature, the final rotational temperature is equal to 0.36. It has to 
be noted that the final distribution is slightly non-Boltzmann (more peaked) with respect to the adjusted thermal 
distribution. Within the present harmonic approximation, it is found that the extent of rotational cooling does not 
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FIG. 6: Distribution of initial and final angular momenta from thermal excitations of LJ„+i at T = 0.5. The best fit to a 
Boltzmann distribution is also reported for the product, (a) n = 50; (b) n = 100. 

strongly depend on cluster size. It is worth comparing the present approach to the results obtained by Stace on Ar+ 
clusters. 4 The physical conditions chosen by this author are those of a cold thermal distribution of J, with an infinitely 
narrow extra vibrational energy distribution located at AE. Because AE is very large in this work (1 eV), many 
evaporations will occur, and the width of the vibrational energy distribution does not really matter. Hence we can 

safely take T v ib = Trot- We have calculated the ratio /T^ot between the final and initial rotational temperatures, 
for LJ50 and at several excitation energies in the range < AE < 80 LJ units. The value used by Stace corresponds 
approximately here to AE ~ 50-80 units. The results are plotted in Fig. [7| At small extra excitation energies, J c is 
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FIG. 7: Final/initial rotational temperature ratio following evaporation of LJ51, for two thermal distributions at T = 0.05 and 
T — 0.20, as a function of the extra excitation energy. 

small and evaporation is more likely accompanied with rotational cooling. On the other hand, at large values of AE, 
rotational heating occurs easily. It seems also pretty obvious that the smaller the initial temperature, the stronger the 
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heating at a given excitation energy. The transition between cooling and heating occurs near AE = 12 (for T = 0.05) 
and AE — 25 (for T = 0.20), and we do not find any evidence for rotational cooling in LJ 50 at excitation energies 
as large as 50 LJ units. This is in clear contradiction with the results by Stace^ who found heating after a single 
evaporation in clusters smaller than 30 atoms, and cooling in 40- and 50-atom clusters at T ~ 0.05. At T ~ 0.20, 
only the 10-atom cluster undergoes rotational heating. This trend is also observed in our work, since an increase in 
the initial rotational temperature indeed favors cooling. 

IV. DISCUSSION AND CONCLUSION 

Systematic investigational have shown that phase space theory, in the sense of Chesnavich and Bowers^ performs 
much better than RRK or Weisskopf-Engelking models in reproducing the data of reference trajectory simulations. 
Weerasinghe and Amar emphasized that the interaction between an atom and the remaining cluster could not be 
modelled by a simple C/r n potential, and that the geometrical extent of the system had to be taken into account. 
The detailed studies by Peslherbe and Hasc 8 on rotating and nonrotating aluminium clusters confirmed the general 
agreement between PST and MD, but they noticed some deviations at large excitation energies. Probably the simple 
form chosen for the dissociation potential is not appropriate for these clusters, for which the interatomic potential 
(Lennard-Jones pairwise + Axilrod- Teller three-body component) is also not fully isotropic. 

To our knowledge, there has not been any attempt at describing evaporation in rotating clusters using PST, except 
in the works by Peslherbe and HaseA& One of the interesting results found in the present work concerns the shape 
of the distribution of the product angular momentum, which was seen to have a discontinuity in its derivative for 
J = J r . The statistics gathered by Peslherbe and Hase [100 trajectories for a given (E, J)S] are probably not sufficient 
to confirm this prediction, and we see that 15 sets of 5000 trajectories are necessary to get a hint at this effect. 

Here, we have used PST in its general form, assuming that both energy and angular momentum were conserved 
during dissociation, calculating transition state properties at the centrifugal barrier, and without assuming a C'/r n 
form for the interaction potential between fragments. The main approximation was to treat the product cluster as a 
sphere, but this seems appropriate, in practice, for a large cluster in its liquidlike phase. 

The significant differences between the present results and the work by Stace^ probably comes from the assumption 
made by this author that all translational energies released are uniformly distributed. This is a very strong assumption, 
which, as a matter of fact, prevents a rigorous conservation of energy during dissociation. Its consequence is that 
the most probable value of St is (E — Eq)/2, much too large with respect to the simple RRKM estimate given by 
[E — Eq)/s. Therefore rotational energies are significantly underestimated, which is precisely observed under the form 
of preferential cooling. To confirm these differences between our work and the results by Stace, we have performed 
additional MD simulations of LJio, at the same conditions as above for n = 50. In particular, we have carefully 
looked at the distribution of final angular momentum of the product LJg, but we could not find any signature of 
bimodality, contrary to the picture discussed in Ref. 0- Our PST calculations, within the harmonic and simple 
potential approximations, show a stronger rotational heating than the Stace prediction: for T r ^ = 0.05 we get 
T'rot = 0.28 at an extra excitation energy AE = 25 units only (in agreement with MD), instead of T^ t « 0.22 for 
AE > 50 units in Stace's worki^ 

Even though the harmonic approximation may be questioned at high energies, it allowed us to undertake quantitative 
studies for clusters containing an arbitrary number of atoms. This required some elementary properties such as 
rotational constants or binding energies to be approximated in a liquid drop fashion, with explicit functions of 
size. Surely this approach is too simplistic, as it neglects non-monotonic finite size effects^ The form chosen for 
the interaction potential may also not be appropriate for very large sizes, for which the location of the centrifugal 
barrier becomes of the same magnitude as the cluster radius itself. Fortunately the numerical effort involved in the 
calculations of the anharmonic DOS and the effective dissociation potential is comparatively smaller than running a 
series of trajectories, especially at low energies where the dissociation rate is vanishingly small. 

The work carried out in this paper and in our previous article^ can be straightforwardly extended to the case of 
multiple (sequential) evaporations, which arise at large excitation energies, and it should be possible to overcome most 
of the limitations that Stace had to imposed Possible connections with the search for a liquid-gas transition in finite 
systematise could then be investigated. Another extension to molecular clusters would also be of great interest, as 
there are presently only very few available theoretical results based on statistical analyses at the level of phase space 
theory. 
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